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Abstract. Within the framework of hierarchical cluster- 
ing we show that a simple Press-Schechter-like approx- 
imation, based on spherical dynamics, provides a good 
estimate of the evolution of the density field in the quasi- 
linear regime up to £ ~ 1. Moreover, it allows one to 
recover the exact series of the cumulants of the probabil- 
ity distribution of the density contrast in the limit £ — > 
which sheds some light on the rigorous result and on "fil- 
tering" . We also obtain similar results for the divergence 
of the velocity held. 

Next, we extend this prescription to the highly non- 
linear regime, using a stable-clustering approximation. 
Then we recover a specific scaling of the counts-in-cells 
which is indeed seen in numerical simulations, over a well- 
defined range. To this order we also introduce an explicit 
treatment of the behaviour of underdensities, which takes 
care of the normalization and is linked to the low-density 
bubbles and the walls one can see in numerical simula- 
tions. We compare this to a 1-dimensional adhesion model, 
and we present the consequences of our prescription for the 
power-law tail and the cutoff of the density distribution. 
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1. Introduction 

In the standard cosmological model gravitational struc- 
tures in the universe have formed by the growth of small 
density fluctuations present in the early universe (Peebles 
1980). These perturbations may be due to quantum me- 
chanical effects and are likely to be gaussian, so that they 
are described by their power-spectrum. In many cases, 
like the CDM model (Peebles 1982; Davis et al.1985), the 
power increases at small scales which leads to a hierarchi- 
cal scenario of structure formation. Small scales collapse 
first, building small virialized objects which merge later to 



form broader and broader halos as larger scales become 
non- linear. These halos will produce galaxies or clusters 
of galaxies, according to the scale and other physical con- 
straints like cooling processes. Hence it is of great interest 
to understand the evolution with time of the density field 
and the mass functions of various astrophysical objects 
it implies. This is a necessary step in order to model for 
instance the formation of galaxies or clusters, which can 
later put constraints on the cosmological parameters using 
the observed luminosity function of galaxies or the QSO 
absorption lines. 

Within this hierarchical framework an analytical 
model for the mass function of collapsed (or just- 
virialized) objects was proposed by Press & Schechtcr 
(1974) (hereafter PS). Numerical simulations (Efstathiou 
et al.1988; Kauffmann & White 1993) have shown this 
mass function is similar to the numerical results, although 
the agreement is improved by a modification of the usual 
density threshold 5 C used in this model (Lacey & Cole 
1994). However, this prescription encounters some serious 
defects, like the cloud-in-cloud problem studied by Bond et 
al.(1991), and a well-known normalization problem since 
one only counts half of the mass of the universe. Moreover, 
underdense regions are not modelled by this approach, and 
are simply taken into account in fine by a global multipli- 
cation of the mass function by a factor 2 which allows to 
get the correct normalization (this multiplicative factor 2 
was recovered more rigorously by Bond et al.1991 for a 
top-hat in k, but this does not extend to more realistic fil- 
ters). Finally, this model only predicts the mass function 
of collapsed objects, while one may also be interested in 
mass condensations of various density, in order to study 
galaxies or even underdense structures ("voids") for in- 
stance. 

A method to describe the density field itself, and then 
derive the multiplicity functions of any objects, is to con- 
sider the counts-in-cells. These are closely related to the 
many-body correlation functions and were studied in de- 
tail by Balian & Schaeffer (1989) in the highly non-linear 
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regime within the framework of the stable clustering pic- 
ture (Peebles 1980), where the correlation functions are 
scale- invariant. This assumption and the scaling it implies 
for the matter distribution are indeed verified by numer- 
ical simulations (Bouchet et al.1991; Colombi et al.1995) 
and observations (Maurogordato et al.1992). Bernardcau 
& Schaeffcr (1991) and Valageas & Schaeffer (1997) (here- 
after VS) described some of the consequences of this model 
for the multiplicity functions of various objects like clus- 
ters or galaxies. Indeed, a great advantage of this approach 
is that one can study many different astrophysical objects 
(and even "voids") from a unique model which is not re- 
stricted to just-collapsed objects. Moreover, it is based on 
a rather general assumption, which docs not depend on 
the exact details of the dynamics, and VS showed that 
the PS prescription can be recovered as a particular case 
among the possible models this scale-invariant approach 
can describe. Note that since the PS approach is unlikely 
to give a very precise description of the clustering process, 
because of its simplicity and the problems it encounters, 
the scale-invariant method offers the advantage to provide 
a simple and natural way to take into account the possible 
corrections to the PS prescription. 

In this article, we intend to show how a simple PS-like 
model, based on the spherical dynamics, can describe the 
evolution of the density field and provide a specific model 
for the scale-invariant picture studied in Balian & Schaef- 
fer (1989) or VS. Thus, wc consider both the quasi-linear 
(E < 1) and highly non-linear (E 3> 1) regimes (where as 
usual E(M, a) is the amplitude of the density fluctuations 
at scale M given by the linear theory at the considered 
time), to relate the final non-linear density field to the 
initial conditions and to show how the scale-invariance of 
the many-body correlation functions and the scalings of 
the multiplicity functions described in VS can arise from 
the hierarchical structure formation scenario. We focus on 
the counts-in-cells, which provide a powerfull description 
of the density field and can be used to obtain any mass 
function of interest as detailed in VS. First, we derive the 
statistics of the counts-in-cells which our spherical model 
implies in the quasi-linear regime (Sect. 2). We show in 
particular that we recover in the limit E — > the whole 
series of the cumulants derived rigorously by Bernardeau 
(1994a) from the exact equations of motion. We also con- 
sider the predictions of this model for the statistics of the 
divergence of the velocity field (Sect. 3) and explain why 
this simple spherical dynamics works so well in the limit 
£ — > (Sect. 4). Finally, we study the non- linear regime 
(Sect. 5). In addition to the virialization process of over- 
densities we take care explicitely of the evolution of un- 
derdensitics. This new prescription solves in a natural way 
the normalization problem (in the sense that all the mass 
will eventually get embedded within overdense halos) and 
is compared with a 1-dimensional adhesion model. It also 
allows to complete the comparison with the scale- invariant 
approach. Indeed, although we do not obtain as detailed 



results for the mass function as the PS approach, which 
we believe anyway to be rather illusory since such simple 
and crude descriptions cannot provide rigorous and exact 
predictions, we think our model allows to understand how 
scaling properties appear, and it can predict asymptotic 
behaviours like the slope of the power-law tail of the mass 
function or its exponential cutoff (so that the model can be 
tested). Moreover, it gives a usefull reference which would 
enable one to evaluate the magnitude of the effects ne- 
glected here from a comparison with the results obtained 
by a more rigorous calculation. 

2. Quasi-linear regime: density contrast 

As VS showed, it is possible to extend the usual PS ap- 
proach to obtain an approximation of the density field and 
the counts-in-cells it implies, in addition to the mass func- 
tion of just-virialized objects one generally considers. Al- 
though VS focused on the highly non-linear regime (using 
the stable-clustering ansatz), we shall here first consider 
the simpler case of the quasi-linear regime E < 1. 

2.1. Critical universe: = 1 
2.1.1. Lagrangian point of view 

We first consider a Lagrangian point of view, which is 
well suited to PS-like approaches since the fundamental 
hypothesis of such approximations is that it is possible 
to follow the evolution of fluid elements recognized in the 
early linear universe (one usually only considers overden- 
sities) up to the non-linear regime. Thus, the usual PS 
prescription assumes i) that the dynamics of these objects 
is given by the spherical model and ii) that their keep their 
identity throughout their evolution. Hence, to any piece of 
matter identified in the early linear universe one can assign 
at a later time a specified radius and density. Of course, 
this cannot be exact as all particles cannot follow simul- 
taneously a spherical dynamics, and this approach cannot 
take into account mergings which change the number of 
objects and destroy their identity. However, we shall see 
below that it can provide a reasonable estimate of the 
early density field. 

In this Lagrangian approach, we consider the evolution 
of "objects" of mass M, identified in the early linear uni- 
verse. In other words, the filtering scale is the mass and 
not the radius. According to the spherical model, parti- 
cles which were initially embedded in a spherical region of 
space characterized by the density contrast Sl (Sl is the 
density contrast given by the linear theory at any con- 
sidered time: Sl oc a where a(t) is the scale factor) find 
themselves in a spherical region of density contrast 5 at 
the same scale M, given at any epoch by: 

5 = ?(5 L ) (1) 
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The function T is denned by the dynamics of the spherical 
model (Peebles 1980): 



1 + 6 = - 



9 (0-sm6) 2 



S L >0 



and 



Sl < : < 



2/3 



1 



Sl 



2 (1- cos 9) 3 



— [6(e -sine)] 



9 (sinhr/ — rj) 2 
2 (cosh 77- l) 3 



"20 [6(sinh?y-?j)] 



(2) 
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This definition of T(8l) breaks down for Sl > S c , with 
(5 C = 3/20 (12tt) 2 / 3 ~ 1.69, where S becomes infinite. 
This is usually cured by a virialization prescription: the 
halo is assumed to virializc in a finite radius taken for in- 
stance as one half its radius of maximum expansion, so 
that F(5 C ) = A c with 1 + A c = 18tt 2 ~ 178. However, we 
shall not consider this modification in this section, as we 
focus here on the quasi-linear regime where most of the 
mass is embedded within regions of space with a density 
contrast S smaller than A c , hence we restrict ourselves to 
S < A c . Now, we define the Lagrangian probability distri- 
bution P m (S): if we choose at random a spherical region of 
mass M its density contrast is between S and S + dS with 
probability P m (6) dS (here the index m does not stand 
for a given mass scale, its use is only to distinguish the 
Lagrangian quantities used here, where we follow matter 
elements, from their Eulerian counterparts which we shall 
introduce in the next sections). Within the framework of 
the spherical model, this probability distribution is simply 
given by: 



P m (6) dS = P L (S L ) dS L with S = T(5 L ) 



(4) 



where Pl(Sl) is the probability distribution relative to 
the linear density contrast, or the early universe. We shall 
assume that the initial density fluctuations are gaussian, 
so that: 



Pl(Sl) = 



1 



-<5|/(2S 2 ) 



2ttE 



(5) 



where as usual E = E(M, a) is the amplitude of the den- 
sity fluctuations at scale M given by the linear theory 
at the considered time (E cx a M~( n+3 ^ 6 for an initial 
power-spectrum which is a power-law: P(k) cx k n ). We 
can note that the probability distribution P m (S) is cor- 
rectly normalized: J P m (5) dS = 1 by definition. However, 
the mean value of 8 is not equal to 0. This is quite natu- 
ral since particles get embedded within increasingly high 
overdensities (while the density contrast cannot be smaller 
than —1), and at late times one expects that all the mass 
will be part of overdense virialized halos. However, this 



late stage is not described by the model introduced in this 
section, which does not include virialization and where 
half of the mass remains at any time within underdensi- 
ties (this is the well-known normalization problem of the 
PS prescription by a factor 2). 

We can also characterize the probability distribution 
P m (S) by the generating function <p m (y,€ m ) (and by £ m ) 
which we define by: 



-<Pm(v,t m )/tr, 



/oo 
-1 



-Sv/Z n 



P m (S) dS 



(6) 



where £ m =< S 2 >. The function ip m (y,^ m ) generates the 
series of the moments of P m (S): 



-<p m (v,i m )/S r . 



P =i 



tll yP < SP > 



Sm 



(7) 



In the highly non-linear regime considered by VS the func- 
tion tpm(y,£, m ) was assumed to be scale- invariant, so that 
it did not depend on £ m , but here this is not the case. If 

< S >~ the usual series of the cumulants is generated 

by </?(?/, ?m) : 

< s >= o : ipmiyru - E y p (8) 



p=2 



pi 



The relation (6) can be written in terms of the linear den- 
sity contrast Sl- 



-<Pm(v,t m )/tr, 
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J — ( 



dS, 



2ttE 
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(9) 



In the quasi- linear regime, that is for £ m — > 0, E — > and 
E 2 /£ TO — > 1, we can use the saddle point method to get: 



■<p m (y) =yT{5 y ) + S 2 y /2 



Sy = 



-yr'(Sy) 



(10) 



where <p m (y) is the limit of <p m (y,£m) f° r Cm ~^ 0- ^ wc 
define T rn = —S y and ^ m (r m ) = T(8 y ), we obtain 



fm(y) =yGm(T m ) - y T m g' m (T m )/2 

r m = -y G' m (T m ) 



(11) 



This is exactly the result derived rigorously by Bernardeau 
(1994a) from the equations of motion, when the matter 
is described by a pressure-less fluid. Hence our approach 
should give a good description of the early stages of grav- 
itational clustering, since it is leads to the right limit for 
the cumulants in the linear regime, which was not obvious 
at first sight. Moreover, it provides some hindsight into 
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the result of the exact calculation, as the function T or 
Q m which represents the spherical dynamics appears nat- 
urally in the expression of ip m (y). We shall come back to 
this point in Sect. 4. We can note that our model implies 
in addition a specific dependence of ip m (y,S, m ) on £m ( or 
S), which could be computed from (9). 

2.1.2. Eulerian point of view 

For practical purposes, one is in fact more interested in the 
Eulerian properties of the density field, where the filtering 
scale is a length scale R. For instance, a convenient way to 
describe the density fluctuations is to consider the counts- 
in-cells: one divides the universe into cells of scale R (and 
volume V) and defines the probability distribution P(5) of 
the density contrast within these cells. Hence, we need to 
relate this Eulerian description to the Lagrangian model 
we developed in the previous section. As was done in VS, 
we shall use: 



(l + S') P(S') dS' = F m (>5,M) 



(12) 



with M = (1 + 8) ~p V . Thus the mass embedded in cells 
of scale R with a density contrast larger than 5 is taken 
equal to the mass formed by particles which are located 
within spherical regions of scale M with a density contrast 
larger than 6. Using our Lagrangian model, we get: 



F m {>5,M) = F mL {>8 L ,M) with 5 = T{5 L ) 



(13) 



where F m L is the mass fraction obtained in the linearly 
extrapolated universe. From (12) and (5) we have: 



(1 + 6) P(5) = - 



8_ 

ds 



-u 2 /2 



dv 



S L /S(R m ) 



2tt 



(14) 



with: 

S = T[5 L ] 
Rt = (1 + 5)R 3 
Finally, we obtain: 



(15) 



P(5) 



1 



1 



8 



2tt 1 + 5 85 



r s L - 






_T,(R m )_ 


exp 


[-mi 



(16) 



as in VS. Of course we recover all the mass of the uni- 
verse J(l + 6) P(8) dS = 1. In fact, half of the mass is 
in ovcrdensities and the other half in underdensities, as 
was the case in the Lagrangian description. However, in 
general this probability distribution is not correctly nor- 
malized: J P(5) dS ^ 1 (but in the linear regime, £ — > 
and £ — > 0, its normalization tends to unity). 



We can still define a generating function tp(y,£), as in 
(6), which leads to: 



/oo 
-1 



dS 



8 



2^ 1 + 5 85 



5 L 



x exp 



1 (z 6 l £ 



(17) 



Using 5 = T(5 L ) = G m (T m ), with r TO = -8 L , we define: 
Ar)=G m [r ^ j 

(18) 

t(S) = 



£[(l + 0(r))V3 R] 



Hence <? m (T m ) = Q(t) = 5. In the quasi-linear regime, 
that is for £ — > 0, £ — > and £ 2 /£ — > 1, the saddle point 
method leads to: 



>(y) =yg(r) - yrg'{T)/2 
t =-vQ'(t) 



(19) 



where ip(y) is the limit of tp(y,£) for £ — > 0. Thus, 
once again we recover the result derived rigorously by 
Bernardeau (1994a). Note that we should have modified 
P(5) in (17) since it is not correctly normalized, however 
if this modification only consists of a multiplication fac- 
tor (which must go to unity in the limit £ — > 0) and a 
change of the shape of P(5) for density contrasts much 
larger than £ as £ — > 0, this does not modify the function 
<p(y) we obtained. If the power-spectrum is a power-law 
P(k) oc k n , (16) can be written: 



P(5) = 



(1 + *) 



(n-3)/6 



27r£(i?) 



x exp 



+ 



6 1 + 5 



(l + S) 



(3+n)/6 



(20) 



If the power-spectrum is not a power-law, our method is 
still valid and we have to use (16). However, for a very 
smooth P(k), like a CDM power-spectrum, an easier and 
still reasonable approximation is to use (20) where n is 
the local slope of the power-spectrum at scale R. We 
shall compare this approximation to numerical results in 
Sect. 2. 3. The fact that we recover the exact generating 
function <fi(y,!;) for £ — > suggests again that our pre- 
scription could provide reasonable results in the early lin- 
ear universe when £ < 1. We shall see below that it is 
indeed the case, by a comparison with numerical simula- 
tions. We can note that our probability distributions (20) 
look rather different from those obtained by Bernardeau 
(1994a) since the high-density cutoff is usually different 
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from a simple exponential. However, they both agree with 
numerical results for 5 <C A c . In fact, neither of these ap- 
proaches should be used for large density contrasts S > A c 
where shell-crossing and virialization play an important 
role. 

2.2. Open universe: < 1, A = 
2.2.1. Lagrangian point of view 

In the case of an open universe, we can still apply the 
method described previously for a critical universe but 
there is now an additional time dependence in the relation 
Sl — S. Thus, (4) becomes: 



P m {5) dS = P L (S L ) d6 L with S = F(5 L ,a) 
In the specific case where A = 0, one defines: 



rjb = cosh 
and 



1 



_ 3sinh?? b (sinh?7 b - r/ b ) _ 
[) ~ (coshr/6-1) 2 



(21) 



(22) 



(23) 



which is the growing mode of the linear approximation 
normalized so that D(t — > oo) = 1. Then, the function 
T(5 L ,a) is given by: 



1 



and 



{ cosh rjb — 1 
V 1 - cos 61 

1 + 



9 -sin 6 
sinh n b - rj b 

■ a \ 2 / 3 



sinh rjb - rjb 



1 + S 



( cosh rjb 
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sinh 77 — 77 

\_ cosh T} — 1 J \ sinh r/b — r/b 

/ • u \ 2/3 

/ sinh 77 — 77 x 



\ sinh 776 - Vb 



In the case f2 — > (rjb — > 00, D(t) — > 1), large over- 
densities have already collapsed, and we are left with: 

/ 2 V 3/2 
Q = : 1 + 5 = h--S L \ (24) 

This simple form for F(8l) leads to: 
(1 + 5)- 5 / 3 



Pm(S) 



2tt£(M) 



exp 



9 

8E2 



(l-(l + <5)- 2 / 3 )' 



(25) 



which provides a convenient estimation for P m (S). As for 
a critical universe we can still define a generating func- 
tion fm(y,^ m ) which allows us to recover the results of 
Bcrnardeau (1994a). 



2.2.2. Eulerian point of view 

Naturally, we can obtain an approximation for the proba- 
bility distribution P(S) of the density contrast within cells 
of scale R in a fashion similar to what we did for a critical 
universe from the Lagrangian probability P m (5). In the 
case of a power-spectrum which is a power-law, and in the 
limit £1 — > 0, we get for instance: 



P(S) 



(1 + ,5)(«- 9 )/ 6 



n + 3+(l-n)(l + 6) 



-2/3 



x exp 



W2itY,(R) 

9(1 + a n )(™+ 3 )/ 3 / , \ ; 



(26) 



As we shall see in the next section, this very simple for- 
mula provides in fact a good approximation to P(S) for all 
values of f2 of interest (even for f2 = 1) due to the weak 
dependence of !F(Sl, fi) on f2. As for a critical universe we 
can also define a generating function <p(y,i;) and recover 
the results of Bernardeau (1994a). 

We can note that Protogeros & Scherrer (1997) ob- 
tained similar results with an approach close to ours. They 
considered "local Lagrangian approximations" where the 
density contrast at the Lagrangian point q, time t, is re- 
lated to its initial value by S = !F(5l). They used sev- 
eral approximations for T including the simplified spher- 
ical collapse model (24) and obtained the probability dis- 
tribution (26), which they modified by introducing an 
ad-hoc multiplicative function N(t) within the relation 
6l — » (1 + 5) in order to normalize properly P(6). How- 
ever, our model differs from theirs by some aspects. Thus, 
our Lagrangian probability distribution is defined from the 
start with respect to a given mass scale M (which might 
be seen as a "smoothing" scale). This appears naturally 
in our approach, and it ensures we always work with well- 
defined quantities. Indeed, for a power-spectrum which is 
a pure power-law with n > —3 the "unsmoothed" den- 
sity field is not a function but a distribution. In fact, one 
cannot characterize a point by a finite density, without 
specifying the scale over which this density is realised. 
Then, the change from the Lagrangian to the Eulerian 
view-point is quite natural, and it provides an Eulerian 
distribution function which differs from the Lagrangian 
one in a very simple and physical manner and which de- 
pends on the power-spectrum. No smoothing procedure 
needs to be applied in fine in order to compare with obser- 
vations: a filtering scale (M or R) is always automatically 
included in our approach. Our approximation also shows 
clearly the dependence on Q of the functions F(Sl) and 
ip(y). Finally, one can note that contrary to Protogeros & 
Scherrer (1997) we did not normalize our probability dis- 
tribution P(S). In fact, we think such a procedure is some- 
what artificial and may lead to an even worse approxima- 
tion. Indeed, if the normalization problem comes mainly 
from a specific interval of the density contrast where our 
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approximation is very bad, a simple normalization pro- 
cedure will not give very accurate results in this interval 
(since the starting values have no relation with the correct 
ones) while it will destroy our predictions in the interval 
where they were fine. Thus, one can fear such "cure" may 
in fact spread errors over all density contrasts. We shall 
come back to this point below, but we can already note 
that for n = 1 the probability distribution (26) cannot be 
meaningfully normalized since / P(S)dS = oo. 

2.3. Comparison with numerical results 

From the results of previous paragraphs, since the func- 
tions ip(y) obtained with our prescription in the limit 
£ — > are exactly those derived by a rigorous calcula- 
tion we can expect that the probability distribution P(S) 
we get should provide a good approximation in the linear 
regime. Thus, Fig.l and Fig. 2 present a comparison of our 
approximation with the results of numerical simulations, 
taken from Bernardeau (1994a) and Bernardeau & Kof- 
man (1995) in the case of a CDM initial power-spectrum in 
a critical universe. We display our predictions for a critical 
universe (relation (20)) and an "empty" universe (relation 
(26)), for a power-spectrum which is a power-law with n 
given by the local slope of the actual power-spectrum. 



log[P(«5)] 
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n = -1.3 
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12 3 4 



8 

Fig. 1. The probability distribution of the density contrast 
P(S). The solid lines present the prediction of our prescrip- 
tion for a critical universe and a power-spectrum which is a 
power- law (relation (16) or (20)), for various £ and n, while 
the dashed-lines show the corresponding curves for an empty 
universe Q, — (relation (26)). The data points are taken 
from Bernardeau (1994a) and Bernardeau & Kofman (1995) 
and correspond to a numerical simulation with a CDM initial 
power-spectrum in a critical universe (thus n is the local slope 
of the power-spectrum). The density fluctuation £ and n were 
measured in the simulation. 

We can see that our approach leads indeed to satis- 
factory results up to £ ~ 1 given its extreme simplicity. 
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Fig. 2. The probability distribution of the density contrast 
P{8) as in Fig.l but for two different values of £. 



As was noticed by Bernardeau (1992), the f2 - dependence 
of the probability distribution is very weak (the dashed 
lines are very close to the solid lines on the figures), which 
means that the simple expression (26) provides a reason- 
able fit for all cosmological models up to £ ~ 1. However, 
as we can see on Fig. 2, our prescription leads to a sharp 
peak for very underdense regions 5 ~ — 1 which increases 
with £ but does not appear in the numerical results. This 
defect is also related to the fact that our probability distri- 
bution is not correctly normalized to unity. This problem 
is due to the expansion of underdense regions, which ac- 
cording to the spherical dynamics grow faster than the 
average expansion of the universe so that in our present 
model these areas occupy after some time a volume which 
is larger than the total volume which is available, which 
leads to a probability which is too large. Indeed, within 
our approach underdense regions can expand without any 
limit while in reality this growth is constrained by the fact 
that on large scales we must recover the average expansion 
cx a 3 . Thus, undcrdcnsities join together after some time 
and their mutual influence alters their dynamics, which 
we did not take into account. 

A simple way to normalize correctly the probability 
distribution P(S) would be to define the latter from the 
generating function <p(y,£) which one would take identi- 
cal to if(y) for any £: this is the method used success- 
fully by Bernardeau (1994a). However, there is a priori no 
fundamental reason why this should be a particularly ef- 
ficient procedure (except from the mere constatation that 
it works). In fact, as we can see on Fig. 2 we can expect 
most of the problem to come from the peak which de- 
velops at low density contrasts, so that one should keep 
the consequences of the evolution of P(S), and <p(y,£), in 
other ranges of S and simply disregard the predictions ob- 
tained in the vicinity of this peak or introduce a specific 
modification for this interval. This is even clearer on Fig. 3 
where we can see that our approximation can still pro- 
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P(S) 




Fig. 3. The probability distribution of the density contrast 
P(5). The solid lines present the prediction of our prescrip- 
tion for a critical universe and a power-spectrum which is a 
power-law (relation (16) or (20)), for various E and n, while 
the dashed-lines show the corresponding curves for an empty 
universe Q = (relation (26)). The data points are taken from 
Protogeros et al.(1997). 



vide reasonable predictions for S > and S ~ 1 although 
it completely fails for underdense regions. Of course the 
agreement with numerical results improves for smaller S, 
as shown on Fig.l and Fig. 2, and becomes excellent in the 
limit E — > 0. We only show on Fig. 3 the largest values of 
S where our approximation still makes some sense, in or- 
der to present its range of validity and to show clearly for 
which values of 5 (underdensities) it breaks down first. As 
Protogeros et al.(1997) noticed, the problem becomes in- 
creasingly severe as n gets larger. However, our predictions 
work better than those used by these authors because we 



did not normalize our probability distribution (see their 
Fig. 2). This was expected since we noticed earlier that 
J P(S)d8 = oo for n = 1, so that it cannot even be nor- 
malized, but this does not prevent our approximation to 
provide very good results for low S, and when £ — > we 
recover the gaussian on any finite interval of the density 
contrast which does not include 5 = — 1. 

3. Quasi-linear regime: velocity divergence 

The prescription we described in the previous paragraphs 
can also predict the statistical properties of the divergence 
of the velocity field. Note that we cannot get any reliable 
result for the shear, since our model is based on a pure 
spherical dynamics, so that the shear is zero for all the 
individual regions of matter we consider. This is clearly 
an important shortcoming of this simple approximation, 
however in the linear regime where density fluctuations are 
small the rotational part of the velocity field decays so that 
after a long time (but still in the linear regime) keeping 
only the growing mode the velocity can be described by its 
divergence, or a velocity potential. Thus, our model does 
not contradict a priori the properties of the linear regime, 
which is of course a first requisite. We define the peculiar 
velocity v by: 



r = u = ax + ax = -r + v 

a 



(27) 



and the divergence 6 (V = d/dx in comoving coordinates, 
and V r = d/dr) by: 



6 



Vv a „ 
— = - V r v 

a a 



(28) 



We can note that in the linear regime, where we keep only 
the growing mode, we have: 



(29) 



where D(t) is the growing mode of the density fluctuations 
and f(Q) ~ n" 6 (see Peebles 1980). 

To use the method we described for the density con- 
trast 6, we must now link the divergence 9 of the velocity 
field to the linear density contrast Sl of our Lagrangian 
elements of matter. One possibility is to make the approx- 
imation that the density is uniform over these individual 
regions, and to use the continuity equation (in physical 
coordinates): 



(30) 



An alternative is to consider the mean divergence over the 
Lagrangian matter element V, which can be expressed in 
terms of the expansion of this global volume (velocity of 
the outer boundary): 



d < 6 > 



i^(V.v)d»* = i/ fl V.dB (31) 
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log[P(0)] 

This last "definition" of is the most natural as it does 
not need any additional information on the density profile. 
However, in both cases we obtain: 



dln{\ + 6) 
din a 



(32) 



which relates 6 to Sl through the relation S = !F(Sl,cl) 
used in the previous sections (note that the derivative in 
(32) is to be understood at fixed v: one follows the motion 
of a given fluid element, so that Sl is also a function of a). 
Then, we obtain the Lagrangian probability distribution 
P m (S) d6 = —Pl(Sl) d$L and its Eulerian counterpart, 
which is simply given by P(6) dO = —P{8) dS (we intro- 
duced a negative sign because d9/d8 < 0). Thus the latter 
can be written: 



P(6) 



-1 1 



d 



2¥ 1 + 5 86 



\ 6l ' 








exp 


[4(1)1 



(33) 



We can also introduce the functions (p m e{y,£, m e) an d 
(pe(y,£,e) (with £ =< 8 2 >). As was the case for the den- 
sity contrast we recover in the limit £ e ~ > the functions 
derived by Bernardeau (1994a). 

If we make the approximation T{8l,o) — Fo(Sl) 
where is the function obtained in the limit CI — > 
(see (24)) we get: 



(34) 



= -f(Q) S L [1- -5 L 



Then, if we define w = 6/f(Q) we obtain: 

P(6) dO = P(w) dm 

with 



(35) 



P{w) 



1 



2ttS 



i 2 

1 m 

3 



(n-ll)/4 



x exp 



~2£2 



, 2 

1 TZ7 

3 



6 

(n-l)/2- 



(36) 



Note that the probability distribution of the reduced vari- 
able w does not depend any longer on the cosmology (we 
only neglected the slight dependence on O of the function 
T). 

We compare this approximation (36) to the results of 
numerical simulations taken from Bernardeau et al.(1997) 
on Fig. 4. We can see that we match the high cutoff of 
P(9), which corresponds simply to the expansion rate of 
a void (zero density), but the negative tail is not well re- 
produced for the case of the open universe. Thus, it seems 
that the description of the velocity field is more sensi- 
tive on the approximations involved in our method than 
the density field. This may be due to the influence of the 
shear, which implies that the velocity can no longer be 



' W i 

Si 1 






\ n = i 


*•// 

/ / 




^ = 0.92- 
n = -\ 


y = 0.2 






- ,/ vC = 0.74 
\y n=-l 


t 





Fig. 4. The probability distribution of the divergence of the 
velocity field P(6). The solid line is the prediction of our pre- 
scription for a critical universe while the dashed-lines corre- 
sponds to an open universe O = 0.2 , A = 0. The data points 
are taken from Bernardeau et al.(1997) for the same condi- 
tions (filled circles for the critical universe and rectangles for 
the open universe). 



determined by a mere scalar (through the potential or 
the divergence). Moreover, we can note that, contrary to 
the Zeldovich approximation for instance, our prescrip- 
tion does not provide the location of particles and their 
velocity field from initial conditions (all regions of space 
cannot follow a spherical dynamics at the same time), it 
only gives an estimate for some probability distributions 
without considering the consistent dynamics of all par- 
ticles simultaneously. Thus, the goal of our approach is 
more modest than such a global modelization and it is 
not entirely self-consistent. However, it appears that after 
accepting this shortcomings we obtain nevertheless rea- 
sonable results for the density field. In fact, as we shall see 
in the next section, we expect large density fluctuations 
to follow a dynamics close to the spherical model while 
the intermediate areas which connect these regions obey a 
complex non-spherical dynamics but as their density con- 
trast is of the order of X they do not play an important 
role for the global shape of P(S) as long as £ < 1. 

4. Spherical density fluctuations 

As we showed in the previous paragraphs, the spheri- 
cal dynamics provides a good approximation to the be- 
haviour of the density field (and a reasonable description 
for the velocity divergence) in the quasi-linear regime, and 
it even gives the exact generating functions ip(y) in the 
limit £ — > 0. This may look surprising, since all regions of 
space cannot follow simultaneously a spherical dynamics 
as we have already noticed, so that one would expect this 
prescription to be always somewhat different from the ex- 
act results. In fact, as we shall see below, this is due to 
the fact that the limit £ — > (or equivalently £ — > 0) 
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constrains the increasingly rare density fluctuations of or- 
der unity to be spherically symmetric. Then a spherical 
dynamics should naturally lead to correct results. Note 
that the mean evolution before virialization of rare large 
density fluctuations was treated by Bernardeau 1994b. 

Thus, let us define Sl (resp. S' L ) as the mean density 
contrast over a sphere of radius R (resp. R') centered on 
a point O (resp. O'), and we note r the vector 00' . Then, 
the conditional probability P(5' l \5l) dS' L to get a density 
contrast S' L knowing that we have a density contrast Sl in 
the first sphere is simply a gaussian: 



where we defined: 



(37) 



f 

Jo 



dk 47rfc 2 W(kR) 2 P(k) 



E" 2 = dk 4irk 2 W{kR)W(kR') P(jfc) 
Jo 



£' = E(i?') 

sin(fcr) 
kr 



hence £ =< S\ > , £' =< S'l > , £" =< S L S' L > 
and 5 L0 = S L 



'L 

E" 2 



T 2 



L ' i ^ "L U L 

E 2 E' 2 - E" 4 



E 2 



and W(x) — 3(sinx — xcosx)/x 3 is the top- hat window 
function (see for instance Bardeen et al.1986 for proper- 
ties of gaussian random fields). Thus, the mean value of 
S' L is 5lo which only depends on the distance r = |r| from 
the point O where the density contrast is constrained to 
be Sl over R, which is obvious from the symmetry of the 
problem. However, the profile of the density fluctuation 
centered on O is usually not spherically symmetric, be- 
cause of the fluctuations of S' L . If we consider now the 
limit E — > (i.e. the normalization of the power-spectrum 
goes to 0) at fixed Sl, the mean value Slo does not change 
but Eo — > 0. Thus, in this limit the profile of the density 
fluctuation centered on O becomes spherically symmetric 
{S' L (r) = 5lo(M)), an d its dynamics is exactly described 
by the spherical model we used in the previous paragraphs. 
Of course, most of the matter (and space) is formed by 
density fluctuations \5l\ ~ E which are not spherically 
symmetric, but these areas correspond by definition to 
density contrasts which tend to as £ — > and they do 
not influence the shape of the functions ip(y) which de- 
pend on finite values of the density contrast. This can be 
seen for instance on (10) which shows that a finite value of 
y corresponds to a finite value of the density contrast 5 y 
given by S y = —y F'{8 y ) while E — > 0. Thus, our approach 
explains the results of the rigorous calculation of the gen- 
erating functions ip(y) in the limit E — > from the exact 
equations of motion, and why they depend simply on the 
spherical dynamics. Moreover, our model is not restricted 
a priori to E « 1 and the comparison with numerical 
simulations shows it gives reasonable results up to E <~ 1. 



5. Non-linear regime 

We showed in the previous paragraphs that a very sim- 
ple model, based on the spherical dynamics, can describe 
the evolution of the density field up to £ ~ 1. It is very 
tempting to try to extend this model into the non-linear 
regime up to £ ~ 200, which is sufficient to describe ap- 
proximately the subsequent highly non-linear regime with 
the help of the stable-clustering ansatz (e.g. Peebles 1980 
for a description of this latter approximation). However, 
this implies that we take into account other processes like 
virialization, which can only be made in a very crude way 
within the framework of the prescription described in the 
previous paragraphs, so that we cannot hope to get ac- 
curate quantitative results. In fact, this highly non-linear 
regime where the probability distribution of the density 
contrast is governed by the properties of virializcd ob- 
jects is probably beyond the reach of rigorous perturba- 
tive methods based on the expansion of the equations of 
motion given by a fluid description with an irrotational ve- 
locity field. Indeed, the fluid approximation itself breaks 
down after shell-crossing and cannot describe collapsing 
halos, so that one has to use the Liouville equation which 
makes the analysis more difficult. Hence it may still be 
worthwile to consider simple models like the one described 
in this article which could give some hindsight into the rel- 
evant processes. 



5.1. Spherical collapse 
5.1.1. Counts- in-cells 

According to the spherical dynamics large overdensities 
decouple slowly from the general expansion of the uni- 
verse, reach a maximum radius R m , turn-around and col- 
lapse to a singularity when their linear density contrast 
Sl is equal to a critical value S c , with 5 C ~ 1.69 for f2 = 1 
(we shall only consider the case of a critical universe in 
the following). However, one usually assumes that such 
an overdensity will eventually virialize (because the tra- 
jectories of particles are not purely radial) into a finite 
radius to form a relaxed halo. Generally, this virialization 
radius R v is taken to be one half of the turn-around ra- 
dius, from arguments based on the virial theorem, so we 
shall note R v = a/2 R m with a ~ 1. If we assume that 
after virialization this halo remains stable, its density con- 
trast will grow as a 3 , while its "linear" density contrast 
increases as a 2 . Thus, we modify the function S = T{Sl) 
which links Sl to 5 so that: 



Sl > S c 



1 + S 



10 

3 a 



61 



(38) 



At the time of collapse where the linear density contrast 
Sl is equal to S c , the actual density contrast of the halo 
is 1 + S = a~ 3 (1 + A c ) ~ a" 3 178. Within the framework 
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of stable clustering, a similar relation can be obtained for 
£ (see VS): 



£> 200 : £(R) 



10 

3ot£ 



(39) 



where Rl is related to R in a fashion similar to R m 

(1 + OR 3 (40) 



Rl 



which comes from the fact that these approximations are 
based on a Lagrangian point of view where one follows 
the evolution of matter elements of constant mass. One 
can expect ^ a since all regions of space do not follow 
the same dynamics, and in fact the previous model only 
applies to overdensities, which form one half of the volume 
and mass in the early linear universe. Then, in a fashion 
similar to what was done in VS, we obtain using (16): 



P(S) = ^ — 



5 + n a 
6 a% 



-7)/6 



x cxp 



2a\ 



fi\ (5+«)/3* 

V 



(41) 



so that the probability distribution of the density contrast 
satisfies the scaling-law: 



P(5) = ^ h(x) with x = - = - 



(42) 



since in the regime we consider here (S > 200, £ > 200) 
we have (1 + £)/£ ~ 5/£. The variable x is related to the 
usual linear parameter v by: 



at 

while the scaling function h(x) is simply: 



x 2 h(x) 



r (5+n)/6 



2tt 



cxp 



,(5+n)/3 



2a\ 



(43) 



(44) 



In fact, this scaling is characteristic of a much wider class 
of models, defined by the scale-invariance of the many- 
body correlation functions, studied in detail by Balian & 
Schaeffer (1989). Thus, the model we described in this ar- 
ticle, based on the spherical dynamics and a strong stable 
clustering assumption, appears as a simple way to esti- 
mate the scaling function h(x) characteristic of the highly 
non-linear regime. Then, all the analysis developed for this 
general class of density fields can be applied to this pecu- 
liar model. As was noticed in VS, the approach presented 
in this article allows one to recover the PS mass function as 
it is based on the same fundamental idea: one follows the 
evolution of individual matter elements from the linear 
regime, described by gaussian probability distributions, 



into the non-linear regime. Indeed, if in a fashion similar 
to PS we identify the fraction F(> M) of matter embed- 
ded within just-virialized objects of mass larger than M 
with the mass contained in cells of scale R with a density 
contrast 6 > A CQ (with (1 + A CQ ) = a~ 3 
and M = (1 + A ca )pV) we obtain: 



a 



178 



F(> M) = F(> A CQ , R) = F L (S C , R m ) 



(45) 



Here Fl is the fraction of matter computed in the linear 
gaussian field, and we used the fact that in our spherical 
model a density contrast 5l over a scale R m in the initial 
gaussian random field is associated to a density contrast 
5 = F{5l) over a scale R in the actual non-linear density 
field as described in (13) and (15), as in the PS prescrip- 
tion. Then, the mass fraction in collapsed objects of mass 
M to M + dM is simply: 

=-^ F > A c a ,R)dR 



d 



(46) 



ORr 



■F L (> S c ,R m )dR 

n 



Since we have Fl(> 5 c ,R m ) = Pl(> $c, Rn 
as we should the PS mass function: 



we recover 



MM) 



dM 
~M 



S_c 

2tt E 



dlnS 



dlnAf 



-<5?/(2£ 2 ) 



dM 
~M 



(47) 



However, (41) is not strictly equivalent to the usual PS 
prescription, since in addition to the spherical model it 
also relies on stable clustering. One should note that we 
did not multiply our probability distribution by the usual 
factor 2, so that we only recover one half of the total 
matter content of the universe since the previous argu- 
ments only apply to initial overdensities. Hence the scaling 
function (44) verifies / xh(x)dx = 1/2 while this integral 
should be normalized to unity. This implies (at least) that 
the approximate h(x) we obtained cannot be used for any 
x. In fact, following the discussion of Sect. 4, we expect the 
spherical dynamics model we used so far to be valid only 
for extreme events \v\ » 1. Hence the approximate P(S) 
and h(x) we got should only apply to x ^> 1. The normal- 
ization problem of the PS mass function is often "cured" 
by an overall multiplication by a mere factor 2, "justified" 
by the excursion set approach in the case of a top-hat in 
k (Cole 1989; Bond et al.1991). However, this result does 
not extend to other window functions (like the top-hat 
in real space used here) with which for large overdensi- 
ties v 3> 1 the PS mass function does not suffer from the 
cloud-in-cloud problem in the sense that the multiplica- 
tive factor needed to correct for double counting goes to 
1 (and not 2) for large masses (Bond et al.1991; Peacock 
& Heavens 1990; VS). Hence, it appears that one should 
not multiply (44) by 2, but merely restrict its application 
to x 3> 1. Moreover, we can note that this normalization 
problem is closely related to the behaviour of underdensi- 
ties, which are not well described in the usual PS approach 
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and constitute the missing half of the matter content of 5.2. Under densities 

the universe. 5.2.1. Density profile of underdensities 



5.1.2. Density profile of virialized halos 

We can notice that an overdensity with an initial den- 
sity profile which is exactly given by the mean value Sl 
(Sect. 4) leads to a final virialized halo with a flat slope in 
its inner parts, since the initial density contrast S' L con- 
verges to a finite value at the center (R' — > in the case 
r = 0). This is not consistent with the results obtained 
from numerical simulations which find inner density pro- 
files p oc (Navarro et al.1996; Navarro et al.1997; Tor- 
men et al.1997) or even steeper (Moore et al.1998). How- 
ever, for moderate values of v the correlation between the 
density at scale R and the density at a smaller scale R' 
enclosed within the former one quickly weakens as the fluc- 
tuation S <~ X' becomes larger than 5lo (£' diverges for 
small R' while Slo remains finite). As a consequence one 
cannot infer the average density profile of virialized halos 
from 5lo- In fact, it seems more reasonable to consider an 
initial density profile of the form: 



(48) 



which leads to a final density profile for the virialized halo: 



(1 + 5) oc p oc R 7 with 7 = 



3(3 + n) 
5 + n 



(49) 



where x (resp. R) is a comoving (resp. physical) coor- 
dinate. The reasoning below (48) is that if we look at a 
virialized halo of mass M, its characteristic density within 
a smaller sphere of mass M' (centered on the peak) will 
be set by the maximum linear density contrast 5' L realised 
over all spheres of mass M' enclosed within the larger 
matter element M. The value of this maximum S' Lmax 
will scale with M' as as soon as £' 3> S, see (37), 
which leads to (48). This picture also assumes that during 
virialization new collapsing shells which may not be cen- 
tered on the density peak will roughly circularize around 
it. We can note that for — 2 < n < — 1 the slope (49) is 
1 < 7 < 1.5 which is close to what is seen in simulations 
(Navarro et al.1997; Moore et al. 1998a). The same reason- 
ing could also be applied to underdensities, which we shall 
use in the next section. However, the previous arguments 
arc quite crude and a much more detailed analysis would 
be required to get a good description of virialized halos. 
Moreover, within the approach described in the previous 
section the shape of the mean density profile has a rather 
weak meaning since the model implies that a lot of sub- 
structure is present within halos, so that a large object 
can be decomposed as a hierarchy of many smaller peaks 
with larger densities. Note that some simulations (Ghigna 
et al.1998) seem indeed to show that many sub- halos can 
survive within larger objects although not to such a large 
extent as in the model (this might be due to finite resolu- 
tion effects). 



The function x 2 h(x) obtained in the previous section (44) 
seems to show that the exponent of the small x power-law 
tail is u) = (5 + n)/6. However, as we argued above we do 
not expect this scaling function to give reasonable results 
for x <C 1, so that we need to get uo from another point 
of view, which considers explicitely low-density regions. 
Indeed, since we expect the spherical dynamics approxi- 
mation to be valid mainly for rare events we shall shift 
from i^3>lto^<C— 1: that is we now study very under- 
dense areas. Moreover, initially non-spherical underdensi- 
ties tend to become increasingly spherical as they expand 
(contrary to the collapse which enhances deviations from 
spherical symmetry), as seen in Bertschingcr (1985), so 
that a spherical dynamics model could give reasonable re- 
sults. We shall assume in this paragraph that the many- 
body correlation functions are scale-invariant, so that the 
density field is described by the scaling function h(x) in 

the non-linear regime £>lfbr(l+(5)>£ M ; which 
includes very underdense and small x regions. Then, if we 
consider a sphere of radius R with a density contrast 5 
over R, the mean density contrast < 1 + 5" >s on its 
outer shell can be shown to be: 



x < 1 : < 1 + 5" > s = 



(! + <*) 



for low density regions ((1 + S) -C £), where 

3(3 + n) 
7 = 7~, 



(50) 



(51) 



is the slope of the two points correlation function in the 
highly non-linear regime and n is the slope of the initial 
power-spectrum P(k) which we assume here to be a power- 
law. This means that the density profile is locally p oc 
R 7"V(i-w)_ Thus, to get u> we may consider the evolution 
of the density profile of a typical spherical underdensity 
using (50). Let us follow an underdensity with an initial 
profile in the early universe \Sl\ oc X(M) oc M~( n+3 ^ 6 
(see (48)). Its dynamics is simply given by the spherical 
model: 



' R = A(coshr) - 1) 

t = i?(sinh?7 — n) 
and 



with A 3 = QMB 2 



Hence we have: 



2/3 



'Soc \S L 



-3/2 



A oc M^^B 2 ' 3 oc |<y-("+5)/("+3) 



(52) 



(53) 



(54) 
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At late times rj ^> 1, and using (3) we obtain: 
' R ~ A cosh r/ ~ A t/B ~ ^("-iJ/WS+n)) 

(1 + 5) ~ 1/ cosh ry - \S L \~ 3 / 2 



(55) 



which means that (1 + S) oc #3(3+n)/(i-n)_ If we identify 
this exponent with 7^/(1 — w) we obtain eventually: 



5 + n 
6 



(56) 



which is also the value one would have inferred from (44). 
We can note that although the initial density is lower 
in the central regions of the perturbation, which expand 
faster, there is no shell-crossing and the density profile is 
given by (55) for n < 1 (the profile of the initial density 
contrast is less steep than Sl oc R~ 2 ). As we could expect 
from the analysis we developed in the linear regime, we 
can note that the case n = 1 leads to some problems since 
we would have to — 1 . Hence it cannot be described by this 
simple model. We shall only consider —3 < n < 1, which 
corresponds to hierarchical clustering and cosmologically 
relevant power-spectra (but one may expect that n — ► 1 
on very large scales R — ► oo). 

5.2.2. Contact of underdensities 

We shall now develop a slightly more detailed description 
of the fate of underdensities, in order to follow the be- 
haviour of the matter (one half of the total mass) which 
was initially embedded in these low-density areas. This 
means that we have to modify the function T{5l) for 
Sl < too. As we noticed in Sect. 2, underdensities grow 
very fast according to the spherical dynamics which leads 
to an approximate probability distribution of the density 
contrast with an ever increasing normalization (instead of 
unity). In fact, the volume formed by any given range of 
underdensities will eventually outgrow the available vol- 
ume of the universe within the approximation used so far. 
Indeed, we can write (14) and (16) as: 



P{5) dS 



with 



1 



1 



2t7 1 + 5 



- 2 / 2 dv 



(57) 



S(i? m ) 



(l + J)(3+n)/6^ 

S(i2) 



The linear variable v defines the underdensities as it is a 
constant of the dynamics. We can see from (57) that a 
logarithmic interval of \v\ of order unity will occupy all 
the volume of the universe when 



1 



^7T l+<5 



1 



(59) 



Naturally, we cannot keep our model unchanged for later 
times, since the mean volume of these underdensities 



should not grow faster than a 3 after this date. In fact, 
in this picture such a range of negative density fluctua- 
tions first expands following the spherical dynamics un- 
til neighbouring underdensities come into contact and fill 
the entire universe. At this time the universe appears to 
be constituted of low density bubbles of size R and den- 
sity p (with p < p). The matter which was "pushed" 
by these "voids" to form the interface between adjacent 
underdensities gets squeezed and reaches high densities. 
Thus, we shall assume the latter build virialized struc- 
tures of size ~ R and density — p (for instance they may 
have a density — 200 times higher than the density p of 
neighbouring areas). These represent the sheets and fila- 
ments one can observe in numerical simulations (e.g. Bond 
et al.1996) which separate low-density regions, while the 
spherical overdensities described in the previous sections 
are the nodes which form at their intersections (note that 
the latter, corresponding to large positive initial density 
fluctuations, appear first). To see more clearly what this 
would imply for the probability distribution of the density 
contrast (and for the mass functions) we shall simply con- 
sider that underdensities defined by their linear parameter 
v, and their scale R, stop expanding when: 



1 + S 



-u 2 /2 



(60) 



with /3 <~ 1 (that is when they fill all of the universe) 
and keep after this date a constant radius R and density 
p. Thus, as time goes on the universe gets filled with in- 
creasingly large and underdense "bubbles" while a grow- 
ing fraction of the matter progressively forms virialized 
structures (with a characteristic density which becomes 
vanishingly small as compared to the mean density) . Note 
that within such a picture all the mass will eventually be 
embedded in virialized high density objects, so that the 
usual normalization problem is solved in a natural way. 
Thus, we must now use (60) to obtain a relation v — S 
in order to get the probability distribution of the den- 
sity contrast from (57). A negative density fluctuation v 
"stops" when it reaches a density contrast S s on a scale 
R such that £ = S s , defined by (58) and (60). At late 
times, rj 3> 1, the spherical dynamics (see (3)) leads to 
the approximate relation: 



(58) (1 + S) 



20 
27 



is\Y,(R) 



6/(n-l) 



(61) 



hence we obtain: 
f(l + <5 s ) = M e - 2 /2 



27 p(l-n)/6 |j,|(n-7)/6 e (l-n) u 2 /l 



(62) 



Note that the scale R is not the Lagrangian scale R m , 
and vY,{R) ^ Sl- After this "stopping time" t s , the radius 
of the object does not evolve any longer while its density 



P. Valageas: Structure formation: a spherical model for the evolution of the density distribution 



13 



contrast increases as a 3 , hence initial underdensities which 
have already "stopped" verify: 



(! + *) = (! + *.)(,) (fg) 



6/(5H 



(63) 



since E(i?) cx a ( 5+ ™)/ 2 i?-(«+3)/ 2 by definition, where i? is 
the physical radius. This is the relation v — 6 we needed to 
derive the probability distribution of the density contrast. 
We can note that the density contrast is now a function of 
both Sl and S, contrary to the pure spherical dynamics 
case used previously where we had 6 = T(8l)- This is 
in fact a necessary condition to be able to get eventually 
all the mass of the universe within overdense virialized 
structures. Finally, this leads to: 



I 2 V2w 2a 5 



-21nf^ 



-3/2 



(64) 



where uu is defined by (56), x by (42) and we used (39) to 
introduce £. Thus, we recover the scaling-law (42) with the 
same exponent u as previously, with logarithmic correc- 
tions. Such logarithmic terms may indeed exist, but our 
model is probably too crude to give a reliable estimate 
of their importance. As was the case for overdensities, we 
also obtain a relation between the linear and non-linear 
parameters v and x: 



v 2 e~ v / 2 = 



9/3 
2c^ 



(65) 



This scaling (42) of P(8) applies to density contrasts larger 
than the one of underdense regions which are currently on 
the verge of filling the entire universe: (1 + S) > (1 + 6 S ) 
with: 



(1 + S) > (1 + 6,) 
hence 

--1/(1-0,) 

x > £ 



(1^1/(20,-2) I -«/(!-«) 



(66) 



(67) 



which is exactly the limit predicted by a general study of 
the models defined by the scale-invariance of the many- 
body correlation functions. Thus, the approach developed 
in this section "explains" in a natural way both the emer- 
gence of the scaling-law (42) for small x and its range of 
validity. Note that (44) derived from the behaviour of over- 
densities only applied to virialized objects S > 200, and 
was in fact restricted to s > 1 as we argued above. The 
previous considerations also mean that, when seen on co- 
moving scale a: at a time defined by the scale factor a, the 
universe appears to be covered by very underdense regions 
of typical density contrast S s with (omitting logarithmic 
terms): 



as long as one remains in the non-linear regime (£(x, a) 3> 
1 that is (1 + S S ) < 1). The "overdensity" (l + 5 s ) tends to 
for large times or small scales as it should, while the vol- 
ume occupied by most of the matter becomes increasingly 
negligible. 

5.2.3. Adhesion model 

We shall now try to compare the approach developed 
above to a very different point of view: the adhesion model, 
in the peculiar case of 1-dimensional fluctuations (in a 3- 
dimcnsional universe). Indeed, from the picture developed 
in the previous paragraphs we do not expect the adhe- 
sion model to give reliable estimates for the counts-in- 
cells at large values of x > 1, since in this range we should 
count roughly spherical virialized halos which have a ra- 
dius larger than the considered scale R: these correspond 
to the overdensities described by the spherical collapse 
seen in Sect. 5.1. Hence the finite value of the virialization 
radius of these objects plays a crucial role in the final prob- 
ability distribution of the density contrast, which is then 
out of reach of the adhesion model where this radius is 
simply zero. However, this model could give a fairly good 
picture of the filaments and sheets which characterize the 
highly non-linear universe as the transverse thickness of 
these structures, smaller than their length or the radius of 
the neighbouring "bubbles" , does not play an important 
role on the counts-in-cclls realised at these latter scales. 
More precisely, we shall try to evaluate the typical den- 
sity contrast seen on a given scale: this corresponds to the 
maximum of P(6) and to the density contrast S s of the 
"bubbles" which cover the universe. 

We shall first consider 1-dimensional density fluctua- 
tions, since in this case the Zeldovich approximation is 
correct until shell-crossing occurs so that we expect the 
adhesion model to provide reliable results. We still define 
n as the index of the power spectrum < |<5fc| 2 >cx k n so 
that we obtain: 



-1 < n < 3 



£(a;) oc a x^ n+1 ^ 2 



(69) 



This range in n ensures that the density fluctuations in- 
crease at small scales while the fluctuations of the density 
potential grow at large scales, so that we are in the do- 
main of the usual hierarchical clustering. We note x (resp. 
q) the comoving Eulcrian (resp. initial Lagrangian) coor- 
dinate of particles, in this 1-dimensional problem (at early 
times x — ► q) . The relation q — x is given by the adhesion 
model for the displacement field, and the density is simply 
obtained from: 

dq 

dx 



r)(x,a) 



dx 



with i]= (1 + 6) 



(70) 



(l + J 8 )~S(a;,a) 6 /( n - 1 )~a 6 /( n - 1 )a: 3 ( n+3 )/( 1 - n ) (68) v(x,a) = - 



We shall follow some of the notations used by Vergassola 
et al.(1994) and we introduce the reduced velocity v and 
velocity potential ip: 

dip V 



dx 



(71) 



14 



P. Valageas: Structure formation: a spherical model for the evolution of the density distribution 



where V — ax is the peculiar velocity. We shall only con- 
sider the cases — 1 < n < 1 where the initial velocity 
(a — > 0) has homogeneous increments and verifies the 
scale- invariance: 

A > : v (x+\l)-v (x) ^ \ {1 - n)/2 [v (x+l)-v (x)](12) 
while the initial velocity potential satisfies: 



where rj* is the most probable value of 77 (note however 
that by definition the mean value of 77 is simply < 77 >= 1). 
Using (77) we obtain: 



(1 + 6)* (Ax, a) ~ a 



-2/(l-n) 



(Ax) 



(l+n)/(l-n) 



(80) 



A > 



Vo(Ax) ^ A ( 3 -")/ 2 M*) 



(73) 



law 



with Vo(0) = V'o(O) = 0. Here, = means "having the same 
statistical properties" . Then, the overall density over the 
cell [2:1,2:2] is simply: 



We shall now consider the approach based on the 
spherical dynamics we described in the previous sections, 
applied to this 1-dimensional problem in order to com- 
pare its prediction to (80). The equation of motion of a 
1-dimcnsional density fluctuation, of longitudinal size 2R, 
centered on the origin, in an universe invariant by trans- 
verse translations, can be written: 



9 i 2 3 



2 Rb 
" t 2 



(81) 



r]([x 1 ,x 2 ],a) = 



x 2 - x\ 

where the point q(x, a) satisfies: 
2 



[q(x 2 ,a) - q(x!,a)] 



— + aip (q) + xq 



is maximum at q(x, a) 



(74) 



(75) 



where Rb corresponds to the Hubble expansion of the mat- 
ter element: flj, a a and Rj Rb — > 1 as t — > (on the 
other hand the comoving transverse coordinates remain 
constant in time). Thus, x — R/Rb is the comoving co- 
ordinate of the outer front which we normalized so that 
x — > 1 for t — ► 0. We can write (81) as: 



as given by the Hopf-Cole solution of the Burgers equa- 
tion (Hopf(1950), Cole(1951)) obtained from the adhesion 
model. If there is a shock at the Eulerian location x, sev- 
eral Lagrangian coordinates q correspond to the same x 
and we choose the smallest one which we note q m (x,a). 
Then, we have: 



3f £+4ti-2i+2=0 



(82) 



The solution of this equation is simply related to the initial 
conditions by: 



x = q m + av (q m ) 



(76) 



2/3 



and 6l 



2/3 



(83) 



(which is simply the Zeldovich dynamics). Using (73) we 
can check that the density f](Ax, a) over a cell of size A2; 
satisfies: 



where 6l is the linear theory density contrast. Indeed, we 
have (1 + 6) = R b /R = 1/x. Thus, we obtain: (1 + 6) = 
(1 — 5l)~ x , so that at late times: 



v(Ax,a) ^ r/(a- 2 /(" +1 ) Ax, 1) 



(77) 



(1 + <5)«1 : (l + 5)~ \6l\- 1 



(84) 



so that one only needs to consider the time a = 1. The 
highly non-linear regime which is of interest to us here 

corresponds to a — > 00 or A2; — > 0. Numerical simulations (1 + 5 S ) <~ (|f|£(.R)) 2 ^ n_1 ^ 
and theoretical arguments (She et al.1992; Vergassola et 
al.1994) strongly suggest that the Lagrangian map q — ► 
x(q, 1) forms a Devil's staircase and that shock locations 
are dense in Eulerian space (for — 1 < n < 1), which can be 
proved rigorously for n = (Sinai 1992). Hence for almost 
every Eulerian coordinate x a Lagrangian coordinate q m 
exists. Using (76) we can see that for two points 2:1 and 
x 2 at time a = 1 we have: 



The linear parameter v is still given by v = <5i/E(i? m ) 
where now R m = (1 + 6)R so that: 



(85) 



Hence according to the model described in the previous 
section underdensities stop expanding and fill the entire 
universe when they reach the density contrast 6 S such that: 



(86) 



x 2 -xi = q m 2 - q m i + [vo(q m 2) - v (q ml )] 



(78) 
■+ 



In the limits \x 2 — Xi\ — > and \q m2 — q m i\ 
we have from (72) the behaviour \v (q m2 ) — v (q m i)\ ~ 
\q m2 — 5mi|^ 1- ™^ 2 - Since n > — 1 the second term in (78) 
becomes negligible, so that we expect: 



A2; ~ (Aq) 



(l-n)/2 



and V *(Ax,l) ~ (Ax) 2 ^ 1 -^- 1 (79) 



where we omitted logarithmic corrections and we used 
(69). Here, x is the considered comoving scale we noted A2; 
in (80). Thus we recover exactly the behaviour seen above 
within the framework of the adhesion model (80). This is 
due to two effects: i) in such a 1-dimensional problem the 
Zeldovich approximation gives the exact dynamics until a 
shock appears, ii) the formation of these large underden- 
sities corresponds to peaks of the initial velocity potential 
which are global maxima over a large scale but all par- 
ticles located in the final broad low-density areas come 
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from a small Lagrangian region (the peak of ^o) so that 
inner properties are given by local characteristics and the 
spherical model is reasonable (the stop of the expansion at 
S s models the global constraints, related to the fact that 
these peaks are only maxima over a finite scale). In fact, 
both models lead to very similar pictures: most of the uni- 
verse is filled by the expansion of initial low-density peaks 
while most of the mass is squeezed between this under- 
densities in virialized objects (for the spherical model) or 
infinitcsimally thin shocks (in the adhesion model). 

For the usual 3-dimcnsional case, the adhesion model 
does not lead to the same results as the spherical prescrip- 
tion and there would be a break at n = —1 (while it occurs 
at n = 1 for both models in the 1-dimensional case): in- 
deed the previous arguments would give to = (5 + n)/4 
instead of (56) which leads to lower characteristic densi- 
ties (1 + S s ). We think this discrepancy could be due to 
the fact that the Zcldovich approximation does not give 
the exact dynamics any more, even before a shock forms, 
so that the shape and size of structures built at late times 
is not correctly described. Indeed, within this approxima- 
tion the physical coordinate r of a given particle evolves 
as: 



r = a(q + ap) 



(87) 



which means that for a spherical underdensity we obtain 
at late times (1 + 5) <~ t~ 2 while the spherical dynamics 
leads to (1 + S) <~ see (55). Hence the Zcldovich ap- 
proximation overestimates the expansion of voids which 
explains the low (1 + 5 S ) and the high u it gives. Thus, 
the adhesion model appears to confirm our spherical pre- 
scription for 1-dimcnsional fluctuations, where the former 
has a rather firm foundation, which builds confidence in 
our model which predicts moreover scaling-laws which are 
actually seen in numerical simulations. For the generic 3- 
dimensional case, this latter result suggests that our pre- 
scription is still valid, while the adhesion model should 
worsen. 



5.2.4. Non-spherical corrections 

As we noticed above, we expect the spherical dynamics to 
describe extreme events \v\ 3> 1. Moreover, in the same 
way as the "cloud-in-cloud" problem is not very important 
for large overdensities v 1, as noticed in VS, which is 
necessary for the approach developed in paragraph 5.1, 
the corresponding "void-in-void" problem disappears for 
v -C — 1 (since the gaussian density field is symmetric 
under v <-> — v) which allows one to use the prescription 
presented in 5.2.2. However, we can see from (65) that 
large negative values of v correspond to very small x. In- 
deed, with ui ~ 0.5 we can check that v < —4 leads to 
x < 1(T 6 , (l + 5 s ) < 1CT 3 and£ > 10 3 . Thus, the range of 
x which can be studied in current numerical simulations 



0.05 < x < 50 may be too small to recover the power- 
law tail with exponent u) and the exponential cutoff pre- 
dicted by our approach. Hence, non-spherical corrections 
may change the value of uj obtained in numerical simula- 
tions. To get an idea of the magnitude and direction of 
such effects, we can study the case where undcrdensities 
only expand along 1 or 2 directions (planar or cylindrical 
symmetry) while the other axis remain(s) constant in co- 
moving coordinates. Thus the 1-dimcnsional problem con- 
sidered in the previous section corresponds to the expan- 
sion along only one direction, while the spherical model 
represents a growth along all three axis. 

For the 1-dimensional growth, (84) implies that we 
have (1 + S s ) ~ £(M) _1 where we do not consider log- 
arithmic corrections (i.e. \v\ ~ 1). Using R m = R^ and 



E(i2ft) cx ti{Ri 



\(5+n)/6 



we obtain (1 + 8 S ) 



--(5+n)/6 

£ on 



scale Rb (the smallest scale of the underdense regions) . If 
we write this exponent as — w/(l — uj) we get: 



1-D 



5 + n 
ll + n 



(88) 



In a similar fashion, a two-dimensional expansion leads to 
{1 + 6) oc ^-(Vis-i)^ and . 



2-D 



(5 + n)(V13- 1) 



(5 + n)(VT3- 1) + 12 



(89) 



We compare on Table 3 the values obtained in numerical 
simulations to these estimations. Thus, although we re- 
cover the increase of ui with n, non-spherical corrections 
appear to be non-negligible. Hence we expect the value of 
u> measured in numerical simulations, which is presently 
close to the 2-D result, to increase somewhat for smaller x 
at smaller scales where £ is larger, and to get closer to the 
value (56) obtained from the spherical dynamics model. 



Table 1. Exponent ui for various indexes n of the 
power-spectrum. The lines ID, 2D and 3D corresponds to (88), 
(89) and (56). The last two lines present the results of numer- 
ical simulations: Colombi et al.(1997) for A and Munshi et 
al.(1998) for B. 



n 


-2 


-1 





1 


ID 


0.33 


0.4 


0.45 


0.5 


2D 


0.39 


0.46 


0.52 


0.62 


3D 


0.5 


0.66 


0.83 


1 


A 


0.3 


0.5 


0.65 


0.65 


B 


0.33 


0.4 


0.55 


0.7 
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5. 3. General picture 

Thus, in the non-linear regime virialized objects should 
form through two different processes according to our 
model. First, large overdensities with a roughly spherical 
shape collapse as in the PS approach to build high-density 
virialized halos. This corresponds to matter elements de- 
scribed by v ^> 1 and i>1. Second, large underdensities 
expand until they fill the entire universe and see their 
dynamics influenced by their interaction with neighbour- 
ing "voids" . The matter "pushed" by these regions forms 
high-density filaments and sheets at their interface, which 
builds virialized structures of increasingly large scale and 
low density (increasingly lower than the mean density of 
the universe). Obviously the dynamics of the filaments 
and walls cannot be described by a spherical model, but 
our approach takes advantage of the fact that the low- 
density "bubbles" they surround may still be described by 
a spherical dynamics, with the addition of other processes 
to take into account the global constraints which stop their 
expansion. This models matter elements with v <C — 1 
and x -C 1. The behaviour of the intermediate fluid ele- 
ments \v\ <~ 1 is certainly quite complex and depends on 
non-local properties through the influence of neigbhour- 
ing peaks and voids. It is not described by our model and 
corresponds to the transition interval of the scaling func- 
tion h(x) around x ~ 1 from its exponential cutoff to its 
power-law tail. Note that in the regime E » 1, low den- 
sity contrast regions \S\ -C 1 are not described by linear 
theory, even though \S\ is small, because of shell-crossing. 
Thus, our model does not provide as accurate predictions 
as the PS mass function, but it can be tested in numerical 
simulations by studying the density profiles of large halos 
or voids, as well as the asymptotic behaviour of h(x). It 
would also be of interest to follow the evolution of initial 
extreme matter elements \v\ 3> 1. Moreover, we think the 
explicit description of underdensities is a necessary step, 
which was not considered in detail in the PS formulation. 
It ensures in a natural way that all the mass will eventu- 
ally be embedded within high-density virialized structures 
(which occupy a negligible volume), so that there is no 
normalization problem, and it allows to describe the com- 
plex structure formed by low-density bubbles, filaments 
and walls, which is seen in numerical simulations (Cole et 
al.1997; Weinberg et al.1996; Bond et al.1996). One can 
note that according to our model, substructures should 
exist within large overdensities, which is not always seen 
in numerical simulations. However, as Klypin et al.(1998) 
argue this may be due to a lack of numerical resolution, 
moreover substructures do appear through counts-in-cells 
numerical studies. 

In the non-linear regime, we have only considered the 
case of a critical universe so far. However, in a low-density 
universe when O becomes small virialized structures no 
longer form as the linear growth factor D(t) tends to a 



constant: the density perturbations freeze in comoving co- 
ordinate. Thus, on small scales where the density fluctua- 
tions arc large, structures formed early when f2 ~ 1 so that 
the analysis developed above for a critical universe can be 
applied. To get the characteristics of virialized structures 
today at these scales one simply needs to consider these 
early times and then extrapolate until today using the ap- 
proximation (which was used throughout above) that viri- 
alized structures keep a constant scale and density while 
the mean density of the universe decreases as p oc a~ 3 . On 
larger scales, fluctuations will never reach the non-linear 
regime since £ tends to a constant as t — > oo, which is 
nearly reached as soon as O ~ 0.1. Hence, on these scales 
one can simply use the quasi-linear description developed 
in Sect. 2. 2. Thus, one gets a complete picture of the den- 
sity field in a low-density universe (except for a transitory 
range where £ ~ 10) within the framework of the approx- 
imation developed in this article. 

6. Conclusion 

Thus, in this article we have developed a simple model for 
hierarchical clustering based on a spherical dynamics. We 
have shown it provides a reasonable approximation to the 
density field, and the divergence of the velocity field, in 
the quasi- linear regime £ < 1. Moreover, it allows one to 
recover the exact series of the moments of the probabil- 
ity distribution of the density contrast in the limit £ — > 
(contrary to other models like the Zcldovich approxima- 
tion for instance), and sheds some light on the rigorous 
results. Then, we have developed a way to extend this 
model into the non-linear regime, by taking into account 
the virialization of high overdensities (as in the PS ap- 
proach) and also the behaviour of very underdense areas 
(in a way different from the PS prescription). This im- 
plies a particular scaling in x of the counts-in-cells over a 
well-defined range which is indeed verified by the results 
of numerical simulations. Thus, our model deals with the 
evolution of the density field in both linear and non-linear 
regimes, for any cosmological parameters (open or criti- 
cal universe) and provides a simple reference to which one 
could compare the results of a more rigorous treatment. 
Naturally, our approximation should be tested in more de- 
tail with numerical simulations. The density profile of very 
underdense regions, the evolution with time of the scale 
and density of typical "voids" , the asymptotic behaviour 
of the scaling function h(x), should allow one to measure 
the influence of non-spherical corrections and the effects 
of substructure disruption on the density field, which are 
not well described by our model, and precise its accuracy. 
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